{ Titrator 3.0 is a chemical equilibrium calculator intended for use
in aqueous systems with acid/base, complexation, dissociation,
precipitation and electron transfer reactions.

Titrator was designed primarily as a teaching tool for chemistry,
biochemistry, geochemistry and environmental chemistry.  It is purely
an equilibrium calculator, and lacks the database, transport, and
kinetic capabilities of more complex and sophisticated programs.  The
user must provide the stoichiometry and log K value (or enthalpy and
entropy) of each chemical species to be considered.

The algorithm is based upon the MicroQL logic of J. Westall, which in
turn follows the algorithm in the MINEQL program of Westall and FMM
Morel.

Currently documentation is incomplete, but may be useful and includes    
    TiManual        Users manual in .doc and .pdf                     
    TiTutorial      Tutorial based on a calcite system, .doc and .pdf 

Current version 3.0, released under GPL by Steve Cabaniss, Dept. of
Chemistry University of New Mexico, who can be contacted at
cabaniss@unm.edu }

{ Solver -- this unit is a set of mathematical routines which solves
the actual equilibrium problem, as well as associated book-keeping
functions.  The Newton-Raphson algorithm for minimizing the mass
balance error is based upon the MicroQL logic of J. Westall, which
in turn follows the algorithm in the MINEQL program of Westall and
FMM Morel.  The current version works well for purely dissolved
systems, but has problems with systems containing more than one
potential solid phase.}

unit Solver;

interface

uses TitrDefs,             { TitrDefs contains the type and variable }
     Dialogs;              { definitions for the routines given here }

procedure DefaultDefn;     { Initialize the Data Set }
procedure DefaultParams;   { Set temperature, IS and iterative controls }
procedure DefaultMC;       { Set Monte Carlo controls }

procedure SortComp;        { Sorts the components in S by CType }
procedure SortSpec;        { Sort  the species in S by SType }

{ InitArrays moves pertinent data from component and species lists to
  arrays which are used in calculation.  Note that although the
  components and species use log base 10, the computational arrays use
  log base e, requiring a factor of ln(10) for conversion

       Components : Tot     --> Total concentration
                    LogX    --> Initial Guess at log molarity
                    Charge  --> Charge on component

       Species :   LogKZero --> Log K of formation at I = 0
                   A        --> Coefficient matrix
                   SpecChg  --> Charge on species
                   Kadjust  --> (square of SpecChg) - (sum square comp charges)

  NewGuess = True uses the guess values for the initial simulation, whereas
  NewGuess = False keeps the results of the previous calculation}

procedure InitArrays(NewGuess : Boolean);

{ Replace is called after a solution to return calculated values from the
  computational arrays LogX, TotErr and Conc to the component and species lists}

procedure Replace;

{ThermoCorr calculates log K values based on DeltaH and DeltaS of each
reaction at the current temperature.  It returns true only if all species
have deltaH and DeltaS information; otherwise it returns false.
Note that unlike IonCorr, ThermoCorr operates on the log K values in the
component array which is part of the data set, and not on the log K vector.}
function ThermoCorr : Boolean;

{SolidSolve is used to solve equilibrium problems in which one or more solids}
{may precipitate or dissolve. The equilibrium for a given set of solids       }
{is solved using the procedure solve above, and then the solution is checked  }
{versus solid solubility of potentially precipitating or potentially dissolving}
{components. If an existing solid should have dissolved, it is eliminated;    }
{if a solid should have precipitated then it is added to the component list.  }
{This process is repeated until                                               }
{  1) no further modifications are needed or                                  }
{  2) the maximum iterations have been exceeded. }

Function SolidSolve : Boolean;

{********************************************************************}
{   SweepSolve systematicaly varies a single component concentration,     }
{   (stored in the variable Xdata) solving the equilibrium for each value.}
{   It returns True if all problems are solved correctly, False if any    }
{   calculation fails to converge                                         }
{********************************************************************}

function SweepSolve: Boolean;

{********************************************************************}
{   TitrateSolve simulates a titration by solving the equilibrium problem     }
{   for a series of added titrant volumes (stored in the variable AddedVol).  }
{   It returns True if all problems are solved correctly, False if any        }
{   calculation fails to converge                                             }
{********************************************************************}

function TitrateSolve : Boolean;

{***************************************************************************}

implementation


{The procedure DefaultDefn creates a blank data set }
procedure DefaultDefn;
begin
    S.NSpec := 0; 
    S.NComp := 0;
    fillchar(S.Head,sizeof(S.head),0);
    S.Head.Title := 'No data set';
end;


{DefaultParams sets the system parameters to the default values}
Procedure DefaultParams;
begin
   S.FindMu := FALSE;       { user sets ionic strength }
   S.UseMu  := False;
   S.Mu     := 0;           { initialize ionic strength to zero }

   S.ThermoCalc  := False;  { use K values, not delta H and delta S}
   S.Temperature := 298.15; { Temperature is 25 degrees C }

   Toler := 1E-5;           { 0.001% error in mass tolerated }
   ItMax := 100;            { allow 100 iterations before quitting }
end;


{DefaultMC sets the Monte Carlo simulation conditions to the default values}
Procedure DefaultMC;
begin
   EnableMonte := TRUE;  { Allow access to MC routines }
   Monte       := True;  { Perform MC calculations }
   MCReps      := 1000;  { 1000 trials per simulation}
   RSeed       := 1;     { Use 1 as default seed }

   MCSave            := FALSE;
   MonteFileBase     := 'MCResult';
   MCStatSave        := FALSE;
   MonteStatFileBase := 'MCSummary';

   DiffX    := 0.01;      { Difference in log X for derivative calculations}
   SensSave := TRUE;      { Save sensitivity matrix}
   SensFile := 'SensMat'; { Default file name}

   NumBins := MaxBins;
end;


{******************************************************************************}
{*******   Now begin matrix manipulation routines *****************************}

{IntSwap swaps two integer values between variable A and B}
procedure IntSwap(var A, B : Integer);
var i : Integer;
begin 
   i := A;
   A := B;
   B := i
end;


{ DoubleSwap swaps two double precision real numbers }
procedure DoubleSwap(var A, B: Double);
var i : Double;
begin 
   i := A;
   A := B;
   B := i;
end;


{ Routines in this section of the program prepare the data (in the
equilibrium system S) for processing as vectors and arrays and then replace
the results in S after calculations are done.}

{ Sorts the components in S by CType }
procedure SortComp; 
var SortDone : Boolean;    {Is the sort finished?}
   com       : integer;    {Component counter}
   CT        : CompType;   {Component type enumerator} 
   
   procedure CompSwap(i,j : integer);    { Swap two components }
   var dummy : component;   {holds a component during the swap}
      spe    : integer;     {species counter}
      temp   : integer;     {holds stoichiometric coefficient during swap}
   begin
      { Swap component data }
      dummy := S.comps[i];     
      S.comps[i] := S.comps[j];
      S.comps[j] := dummy;
      
      { Now swap columns of coefficients which define species compositions }
      for spe := 1 to S.NSpec do 
      begin
         temp := S.specs[spe].coeff[i];
         S.specs[spe].coeff[i] := S.specs[spe].coeff[j];
         S.specs[spe].coeff[j] := temp
      end;
   end;{ end CompSwap }

begin  { Begin SortComp }
   with S do
   begin
      { Sort with bubblesort }
      repeat
         SortDone := true;
         for com := 1 to NComp -1 do
         begin
            if Comps[com].Ctype > Comps[com+1].CType then
            begin
               Compswap(com,com+1);
               SortDone := False
            end;
         end;
      until SortDone; { End of bubblesort }
      
      {Now count the number of each type of component (first set each
      count to zero)}
      for CT := CTotDiss to CSurface do 
      begin
         CTypeNum[CT] := 0;
      end;
      
      {Increment the # of each component type}
      for com := 1 to NComp do 
      begin
         CTypeNum[comps[com].CType] := CTypeNum[comps[com].Ctype] + 1;
      end;
   end; { end of with S }
end; { end of SortComp }


{ SortSpec sorts the species by class using a bubble sort }
procedure SortSpec;       { Sort species by class }
var SortDone : Boolean;   { Is the sort finished? }
   spe       : Integer;   { Species counter }
   ST        : SpecType;  { Species type enumerator } 
   
   { Swap two species }
   procedure Swap(i,j : Integer); 
   var dummy : specie;
   begin     
      dummy := S.specs[i];
      S.specs[i] := S.specs[j];
      S.specs[j] := dummy
   end;

begin { Begin Sortspec }
   with S do
   begin
      { Sort with bubblesort }
      repeat 
         SortDone := true;
         for spe := 1 to NSpec - 1 do
         begin
            if specs[spe].SType > specs[spe+1].SType then
            begin
               swap(spe,spe+1);
               SortDone := False
            end;
         end;
      until SortDone;

      { Now count the number of each type of species (first set each
      count to zero)}
      for ST := Diss to Sorbed do 
      begin
         STypeNum[ST] := 0;
      end;
      
      for spe := 1 to NSpec do  {Increment the # of each species type}
      begin
         STypeNum[specs[spe].SType] := STypeNum[specs[spe].Stype] + 1;
      end;

   end;
end; { End of SortSpec }


{ InitArrays moves pertinent data from component and species lists to arrays
  which are used in calculation.  Note that although the components and species
  use log base 10, the computational arrays use log base e, requiring
  a factor of ln(10) for conversion

       Components : Tot     --> Total concentration
                    LogX    --> Initial Guess at log molarity
                    Charge  --> Charge on component

       Species :   LogKZero --> Log K of formation at I = 0
                   A        --> Coefficient matrix
                   SpecChg  --> Charge on species
                   Kadjust  --> (square of SpecChg) - (sum square comp charges)}

procedure InitArrays(NewGuess : Boolean); { initialize spec and comp arrays.}
{ If NewGuess is True, put Guess in LogX }
{ If NewGuess is False: keep LogX results }
var  sum, sumsquare : integer;
   com, spe         : integer;            {counters for components and species }

begin
   with S do
   begin
      { fill Tot, LogX, Charge arrays }
      for com := 1 to NComp do 
      begin
         Tot[com] := comps[com].TotalX;
         if (NewGuess or (Comps[com].CType = CFixDiss)) then
         begin
            LogX[com] := Comps[com].GuessX * 2.3025851;
         end;
         Charge[com] := Comps[com].Charge
      end;

      { fill LogK, A, SpecChg, KAdjust }
      for spe := 1 to NSpec do   
      begin
         LogKZero[spe] := specs[spe].logK * 2.3025851;
         sum := 0;
         sumsquare := 0;
         for com := 1 to NComp do
         begin
            A[com,spe] := specs[spe].coeff[com];
            if A[com,spe] <> 0 then
            begin
               sum := sum + round(A[com,spe] * Charge[com]);
               sumsquare := sumsquare + 
               round(abs(A[com,spe]) * Charge[com] * Charge[com])
            end;
         end;
         SpecChg[spe] := sum;
         Kadjust[spe] := sum * sum - sumsquare;
      end;
   end; { end with S }
end; { end InitArrays }


{ Replace is called after a solution to return calculated values from the
computational arrays LogX, TotErr and Conc to the component and species
lists }
procedure Replace;       { return results to comps[] and specs[]}
var spe, com : integer;  { counters for components and species}
begin        
   with S do
   begin              
      {For each component in the system S}
      for com := 1 to NComp do
      begin
         {return the equilibrium }
         comps[com].freeX := logX[com] / 2.3025851;   

         {log free activity}
         comps[com].error := TotErr[com];      {the mass error }
         comps[com].Conc := exp(logX[com]);    {and the free molarity}
      end;
      
      { For each species in the system S}
      for spe := 1 to NSpec do 
      begin 
         {return equilibrium molar concentration}
         specs[spe].conc := conc[spe];     
         if (conc[spe] > 0) then
            specs[spe].LogC := ln(conc[spe])/ LogConvr;
      else 
         Specs[spe].LogC := -999;
      end;
   end;
end; { end of procedure Replace }


{ MatSolve solves the simultaneous linear equations A X = B, where A is the 
coefficient matrix, B is the constant vector and X is the solution vector 
the matrix A is n x n, and the vectors are n units each 
This routine is from Pizer and Wallace, 1983, Chapter 6, and uses  
Gaussian substitution with scaling and pivoting to reduce error }
function MatSolve(var A: matrix; var X,B: vector; n: Integer): Boolean;
var  scale  : vector;  { scaling vector }
   p        : array[1..MaxComps] of integer; { permutation index }
   max      : double;                        { largest scaled pivot }
   maxrow   : integer; { row containing max}
   scaleval : double;  { magnitude of pivots }
   amult    : double;  { row multiplier }
   temp     : integer; { p holder }
   i,j,k    : integer; { counters }
            
begin 
   MatSolve := FALSE;
   
   { Initialize }
   for i := 1 to n do
   begin
      p[i] := i;       { initialize all pivot values}
      scale[i] := 0;   { initialize all scaling values}
      for j := 1 to n do
      begin
         scale[i] := scale[i] + abs(a[i,j]);
      end;
   end;

   { Triangularize w/ pivoting }
   for i := 1 to n-1 do 
   begin
      max := 0;

      { Find the pivot }
      for j := i to n do
      begin
         ScaleVal := abs(A[ p[j],i] / scale[p[j]]);
         if ( ScaleVal > Max) then
         begin
            max := scaleval;
            maxrow := j
         end;
      end;

      { Switch rows }
      temp := p[maxrow];
      p[maxrow] := p[i];
      p[i] := temp;

      { Eliminate subdiagonal elements on the ith column }
      for j := i+1 to n do 
      begin
         if (A[p[i],i] = 0) then 
            Exit;

         amult := A[p[j],i] / A[p[i],i];
         for k := i+1 to n do
            A[p[j],k] := A[p[j],k] - amult * A[p[i],k];

         b[p[j]] := b[p[j]] - amult * b[p[i]]
      end;
   end; { End triangularization w/ pivoting }

   { Back substitution }
   for i := n downto 1 do 
   begin
      x[i] := b[p[i]];
      for j := i+1 to n do
         x[i] := x[i] - A[p[i],j] * x[j];

      if (A[p[i],i] = 0) then 
         Exit;

      x[i] := x[i] / A[p[i],i]
   end; { Back substitution done }

   MatSolve := True;

end; { end Matsolve }




{ MuCalc calculates the ionic strength from species concentration and ionic
charge data }
procedure MuCalc(    NumComp, NumSpec : Integer;
                 var Charge           : IntVector;
                 var FreeX            : Vector;
                 var SpecCharge       : SpecInts;
                 var Conc             : SpecArray);

var   sum : Double;
   i      : Integer;
          
begin
   sum := 0;

   { Calculate the sum of the concentration of each species multiplied by its
   squared charge }
   for i := 1 to NumComp do                   
   begin
      sum := sum + sqr(charge[i]) * FreeX[i]; 
   end;
   for i := 1 to NumSpec do
      sum := sum + sqr(Specchg[i]) * Conc[i];

   { If poor guesses lead to absurd concentrations, set I=1 as max }
   If (sum > 2) then sum := 2;   

   S.Mu := sum / 2;
end; { End of procedure MuCalc }




{ IonCorr corrects the log K values for ionic strength by either the Davies
  equation (I > 0.1) or the Guntelberg approximation (I <= 0.1) }
procedure IonCorr(    NumSpec : integer;          {Number of species}
                  var LogK,                   {calculated Log K at IS = mu}
                      LogKZero,               {given Log K at IS = 0}
                      Kadjust : Specarray);   {Adjustment factor }
var   mult    : double;       { multiplier }
   Davies     : double;
   Guntelberg : double;
   TempFac    : double;          { Temperature dependent factor A }
   spe        : integer;         { species counter }
begin
   {First check to see if IS correction is needed and if IS is high enough to }
   {be important.  If not- simply copy over the existing log K values into }
   {the computational array LogK. }
   
   if ((Not S.UseMu) or (S.Mu < 0.00001)) then
   begin
      LogK := LogKZero;
      Exit;
   end;

   mult := sqrt(S.mu);

   { Temperature factor is proportional to (Kr * T) to the -3/2 power }
   { Kr (relative permitivity is 78.54 }
   { proportionality factor is 1.825 e6 }
   Tempfac := 78.54 * S.Temperature;
   TempFac := 1.825E+6 / sqrt(tempfac * tempfac * tempfac);

   Davies := logConvr * TempFac * (mult / (1 + mult) - 0.2 * S.mu);
   Guntelberg := logConvr * Tempfac * mult / (1 + mult);

   { Guntelberg correction}
   if (S.Mu <= 0.1) then 
   begin
      for spe := 1 to NumSpec do
      begin
         LogK[spe] := LogKZero[spe] + Guntelberg * Kadjust[spe];
      end;
   end
   { Davies }
   else 
   begin 
      for spe := 1 to NumSpec do
      begin
         LogK[spe] := LogKZero[spe] + Davies * Kadjust[spe];
      end;
   end;
end; { end IonCorr }


{ThermoCorr calculates log K values based on DeltaH and DeltaS of each reaction}
{at the current temperature.  It returns true only if all species have deltaH }
{and DeltaS information; otherwise it returns false}
{Note that unlike IonCorr, ThermoCorr operates on the log K values in the }
{component array which is part of the data set, and not on the log K vector}
function ThermoCorr: Boolean;
var spe : Integer;
begin   
   ThermoCorr := True;
   for spe := 1 to S.NSpec do 
   begin
      with S.Specs[spe] do 
      begin
         if KCalc then 
         begin
            DeltaG := DeltaH - S.Temperature * DeltaS/1000;
            LogK := -DeltaG * 1000 /(GasR * S.Temperature * LogConvr);
         end
         else
            ThermoCorr := False;
      end; {end with}
   end; {end for}
end; {end ThermoCorr}


{ The procedure Solve is the main routine for solving an equilibrium problem
  within given mass balance constraints for a given set of components.  The
  species concentrations are calculated from free component activities and
  log K values, then the jacobian matrix Jac is calculated from these
  concentrations and the stoichiometric coefficients in A.  New estimates of
  component activities (G) are then calculated by    solving the system Jac x
  G = TE where TE is the total mass balance vector.  The process repeats
  until either a) the mass balance errors are small or b) the maximum number
  of iterations is exceeded.

  See the MINEQL users guide by Westall and Morel for further illumination. }

{ Params:
    TotalKnown - the # of components for which the total dissolved amount
                 is known.
    NumSpec    - the number of species }
function Solve(TotalKnown, NumSpec: Integer): Boolean;             
var
   {Keeps track of magnitude of summed contributions to mole balance}
   AbsSum   : Double;    
   com, spe : Integer;   {counters for components and species}
   i        : Integer;   {array indices}
   dummy    : Double;    {temporary sum}
   Jac      : matrix;    {Jacobian matrix}
   FreeX    : vector;    {vector of free component concentrations}

begin
   PhaseViol := False;  {No phase rule violation}
   TooMany   := False;  {Not too many iterations}
   Solve     := False;  {not solved}
   MatError  := False;  {Matrix algebra failure flag}
   
   {Begin by correcting the log K values for ionic strength}
   IonCorr(NumSpec, LogK, LogKZero, Kadjust );
   
   Iter := 0; 

   repeat  { Main iterative loop for Solve}

      { Calculate log concentration of each species from the sum of the
      log(K) of formation plus the sum of log(C) of each component times
      the stoichiometric coefficient }
      for spe := 1 to NumSpec do        
      begin                             
         dummy := logK[spe];           
         for com := 1 to S.NComp do    
         begin
            if (A[com,spe] <> 0) then 
            begin
               dummy := dummy + A[com,spe] * logX[com];
            end;
         end;

         { Convert log [] to molarity }
         conc[spe] := exp(dummy); 
      end;

      { Calculate concentration of each component }
      for com := 1 to S.NComp do       
      begin
         FreeX[com] := exp(logX[com]); 
      end;


      { Ionic Strength Calculations }
      if (S.FindMu) then
      begin
         { First find mu based on current concentrations }
         MuCalc(S.NComp, NumSpec, Charge, FreeX, SpecChg, Conc);

         { Then adjust logK accordingly }
         IonCorr( NumSpec, LogK, LogKZero, Kadjust);
      end;

      { Mass Balance (actually mole balance)}
      MassOK := True;

      { Calculate error in total moles of each component}
      for com := 1 to S.NComp do 
      begin 
         {Start with the difference between the free and total for a comp}
         dummy := FreeX[com] - Tot[com];

         {add a summed term to indicate magnitude of total moles}
         {Note: This is needed components with both - and + stoichiometries}
         AbsSum := Abs(Tot[com]) + Abs(FreeX[com]);

         {now add in each species concentration weighted by stoichiometry}
         for spe := 1 to NumSpec do
         begin
            if (A[com,spe] <> 0) then
            begin
               dummy := dummy + A[com,spe] * Conc[spe];
               AbsSum := AbsSum + abs(A[com,spe] * Conc[spe])
            end;
         end;

         TotErr[com] := dummy; 

         { Test relative error }
         if (com <= TotalKnown) then 
         begin
            { Flag if error is too high }
            if (abs(dummy / AbsSum) > toler) then
            begin
               MassOK := False; 
            end;
         end;
      end;

      if (TotalKnown > 0) then 
      begin

         { Compute Jacobian = TotalKnown * TotalKnown }
         for com := 1 to TotalKnown do 
         begin
            for i := 1 to TotalKnown do  { each element dTot / dFreeX }
            begin
               dummy := 0;
               for spe := 1 to NumSpec do
               begin
                  if (A[com,spe] <> 0) and (A[i,spe] <> 0) then
                  begin
                     dummy := dummy + A[com,spe] * A[i,spe] * Conc[spe] / FreeX[i];
                  end;
               end;

               if (com = i) then 
               begin
                  dummy := dummy + 1;
               end;

               jac[com,i] := dummy; { diagonal }
            end;
         end;


         { Calculate new FreeX, LogX }
         { First consider case with only one unknown free component molarity}
         if (TotalKnown = 1) then 
         begin
            NewX[1] := TotErr[1] / Jac[1,1];
            Done := True;
         end
         { If more than one, then use Matsolve to find -deltaX}
         else 
            Done := MatSolve(Jac, NewX, TotErr, TotalKnown);

         if (not Done) then 
         begin
            MatError := True;
            Exit;
         end;

         { Calculate new free X }
         for com := 1 to TotalKnown do  
         begin
            NewX[com] := FreeX[com] - NewX[com];

            { Avoid negative convergence }
            if NewX[com] <= 0 then
            begin
               NewX[com] := FreeX[com] / 10;
            end;

            { Find log(X) }
            LogX[com] := ln(NewX[com]);
         end;
      end; { end TotalKnown if }

      Iter := Iter + 1;  

   until (MassOK) or (Iter = Itmax) or (not Done); { iterative loop complete }

   { Solution converged }
   if MassOK then 
      Solve := True;  

   { Iterations exceeded }
   if (Iter = ItMax) then 
      TooMany := True;

end; {End of Procedure Solve}


{ SolidSolve is used to solve equilibrium problems in which one or more solids
may precipitate or dissolve. The equilibrium for a given set of solids       
is solved using the procedure solve above, and then the solution is checked  
versus solid solubility of potentially precipitating or potentially dissolving
components. If an existing solid should have dissolved, it is eliminated;    
if a solid should have precipitated then it is added to the component list.
This process is repeated until:
  1) no further modifications are needed or                                  
  2) the maximum iterations have been exceeded.} 

{Needless to say, this is a highly empirical and somewhat klugey routine.}
{It seems to work for simple systems, but frequently crashes for more complex}
{systems with many possible solids.}
Function SolidSolve : Boolean;
Const  MaxCount = 10;
Var SolidTries    : integer;  { Number tries at a solution }
   PrecipCount    : integer;  { number of solids precipitated }
   NowSpe1        : integer;  { number of species in current system }
   NowComTot      : integer;  { number of known-total comps after transform }
   NowComDiss     : integer;  { number of known-concentration comps }
   SpeCtr, ComCtr : integer;  { Spec and comp to precipitate/dissolve }
   i              : integer;

SpecSwaps : array[1..MaxSpecs] of integer;   { positions of species and }
   CompSwaps : array[1..MaxComps] of integer;   { components after swapping }


{ TransposeSpec swaps the positions of two species in the arrays LogKZero and
A }

   Procedure TransposeSpec(SwapSpec1, SwapSpec2 : integer);
   Var
      ComCnt : integer;
   begin      
      IntSwap(SpecSwaps[SwapSpec1],SpecSwaps[SwapSpec2]);

      { First, Swap the log K values }

      DoubleSwap(LogKZero[SwapSpec1],LogKZero[SwapSpec2]);
      DoubleSwap(Kadjust[SwapSpec1],Kadjust[SwapSpec2]);
      IntSwap(SpecChg[SwapSpec1],SpecChg[SwapSpec2]);

      { Next, Swap the stoichiometric coefficients A }

      For ComCnt := 1 to S.NComp do
         DoubleSwap(A[ComCnt, SwapSpec1], A[ComCnt,SwapSpec2])

   end;

{ TransposeComp swaps the positions of two components in the arrays Tot and
A }

   Procedure TransposeComp(SwapComp1, SwapComp2 : integer);
   Var
      SpeCnt : integer;

   begin
      IntSwap(CompSwaps[SwapComp1], CompSwaps[SwapComp2]);

      { First, Swap the Tot values }

      DoubleSwap(Tot[SwapComp1],Tot[SwapComp2]);
      DoubleSwap(LogX[SwapComp1],LogX[SwapComp2]);
      IntSwap(Charge[SwapComp1],Charge[SwapComp2]);

      { Next, Swap the stoichiometric coefficients A }

      For SpeCnt := 1 to S.NSpec do
         DoubleSwap(A[SwapComp1, SpeCnt], A[SwapComp2, SpeCnt]);
   end;

{ Transform performs a transformation of basis using component CompTr and
{ species SpecTr. The arrays LogKZero, A and Tot are transformed. }
{ In effect, this routine allows SolidSolve to change from one set }
{ of components to another. }

   Procedure TransformBasis(SpecTr, CompTr : integer);

   Var Pivot       : double;    { coeff of component CompTr in species SpecTr }
      SpecCoeff    : SpecArray; { coeff list for component CompTr }
      CompCoeff    : Vector;    { coeff list for species SpecTr }
      i,j          : integer;   { counters }
      OldK, OldTot : double;    { hold the old values for transformed terms }
      IntTemp      : integer;
      SumSquare    : double;

   begin
      IntTemp := -1 * CompSwaps[CompTr];
      CompSwaps[CompTr] := -1 * SpecSwaps[Spectr];
      SpecSwaps[SpecTr] := IntTemp;

      IntSwap(Charge[CompTr],SpecChg[SpecTr]);

      { calculate Kadjust here }

      { First, save critical values of A for computing new terms }
      Pivot := A[CompTr, SpecTr]; 
      If (Pivot = 0) then begin   
         ShowMessage('Pivot is zero');
         {Need to add error handling code here !!!}
         Exit;
      end;

      OldK := LogKZero[SpecTr];     { Record Log K of species SpecTr }
      OldTot := Tot[CompTr];        { Record Total molarity of CompTr }

      SpecCoeff := A[CompTr];       { Record the species coefficients }

      For i := 1 to S.NComp do      { Record the component coefficients }
         CompCoeff[i] := A[i, SpecTr];

      { Now, perform the transformation of basis by calculating the new values
      of K, Tot and A. }

      LogKZero[SpecTr] := 0;
      For i := 1 to S.NSpec do                               { Transform Log K }
         LogKZero[i] := LogKZero[i] - (SpecCoeff[i] / Pivot) * OldK;

      For i := 1 to S.NComp do                               { Transform Total }
         Tot[i] := Tot[i] - (CompCoeff[i] / Pivot) * OldTot;
      Tot[CompTr] := OldTot / Pivot;

      For i := 1 to S.NComp do
         A[i,Spectr] := 0;
      A[Comptr,SpecTR]  := 1;

      For i := 1 to S.NComp do                               { Transform A }
         For j := 1 to S.NSpec do
            A[i,j] := A[i,j] - (CompCoeff[i] * SpecCoeff[j] ) / Pivot;

      For j := 1 to S.NSpec do
         A[CompTr,j] := SpecCoeff[j] / pivot;

      sumsquare := 0;
      for i := 1 to S.NComp do
         sumsquare := sumsquare + round(abs(A[i,SpecTr]) * Charge[i] * Charge[i]);
      Kadjust[specTr] := SpecChg[SpecTr] * SpecChg[SpecTr] - sumsquare;


   end;          { End of TransformBasis }



{The function DissolutionCheck examines precipitated solids to see if they
{ should dissolve.  The function returns zero if no solid should have }
{dissolved.  Otherwise it returns the number of the dissolved solid}
   Function DissolutionCheck : integer;
   var i, count : integer;            {loop counters}
      Least : double;                {holds most negative mole balance number}

   begin
      Count := 0;
      Least := 0;
      i := NowComTot + 1;
      while ( i <= (NowComTot + PrecipCount)) do begin
         if (TotErr[i] < Least) then begin    { find solid with most negative }
            Count := i;                       { total concentration }
            Least := TotErr[i];
         end;
         i := i + 1;
      end;                                    { End while loop }
      DissolutionCheck := count;
   end; { End of DissolutionCheck }





{The function PrecipitationCheck examines potential solids to see if they }
{should precipitate.  The function returns zero if no solid should have }
{precipitated; otherwise it returns the number of the precipitating solid}
   Function PrecipitationCheck : integer;
   var i, j,  count : integer;   { loop counters }
      Sum, Most : double;     

   begin
      Count := 0;
      Most := 0;
      i := NowSpe1 + 1;
      while ( i <= S.NSpec) do begin    { While loop looks for Ksp most }
         Sum := LogKZero[i];            { exceeded by component }
         for j := 1 to S.NComp do       { concentrations, returns number }
            If (A[j,i] <> 0) then
               Sum := sum + A[j,i] * LogX[j];
         if (Sum > Most) then begin
            Count := i;
            Most := Sum;
         end;
         i := i + 1;
      end;                              { End while loop }
      PrecipitationCheck := count;
   end; { End of PrecipitationCheck}







{ Precipitate causes a solid phase to form via a transformation of basis }
   Procedure Precipitate(Var SpeCtr : integer);
   begin
      { This message was not strictly necessary, but could be implemented}
      {  TopMessage('Precipitation occuring...'); }

      If (SpeCtr <> (NowSpe1 + 1)) then begin       { swap species }
         TransposeSpec(SpeCtr,NowSpe1 + 1);
         SpeCtr := NowSpe1 + 1;
      end;

      ComCtr := 0;
      repeat                                        { find comp in solid }
         ComCtr := ComCtr + 1;
      until ((A[ComCtr, SpeCtr] > 0) or (ComCtr = NowComTot));

      { if comp not found in mass-limited solids, too many solids }

      If ((A[ComCtr,SpeCtr] = 0) and (ComCtr = NowComTot)) then begin
         {This error message really must be implemented}
         ShowMessage('Gibbs Phase Rule violated by solids');
         PhaseViol := True;
         Exit;
      end;

      { Now swap components, if necessary }

      If (ComCtr <> NowComTot) then begin
         TransposeComp(ComCtr,NowComTot);
         ComCtr := NowComTot;
      end;

      { Transform basis by swapping component ComCtr and species SpeCtr }

      TransFormBasis(SpeCtr,ComCtr);

      LogX[ComCtr] := 0;
      PrecipCount := PrecipCount + 1;    { Precipitation }
      NowSpe1 := NowSpe1 + 1;
      NowComTot := NowComTot - 1;
      NowComDiss := NowComDiss + 1;

   end;                    {End of Procedure Precipitate}






{ Dissolve causes a solid phase to dissolve via a transformation of basis }
   Procedure Dissolve(ComCtr : integer);
   begin
      { Another unnecessary message}
      {  TopMessage('Dissolution occuring...');  }

      If (ComCtr <> (NowComTot)) then begin       { swap components }
         TransposeComp(ComCtr,NowComTot);
         ComCtr := NowComTot;
      end;

      {Note- is this right?}
      SpeCtr := S.STypeNum[Diss];


      repeat                                        { find comp in solid }
         SpeCtr := SpeCtr + 1;
      until ((A[ComCtr,SpeCtr] > 0) or (SpeCtr = NowSpe1));

      { if species containing comp not found, error has occurred }

      If ((A[ComCtr,SpeCtr] = 0) and (ComCtr = NowComTot)) then begin
         ShowMessage('Error in procedure Dissolve');  
         PhaseViol := True;
         Exit;
      end;

      { Now swap components, if necessary }

      If (SpeCtr <> NowSpe1) then begin
         TransposeSpec(SpeCtr,NowSpe1);
         SpeCtr := NowSpe1;
      end;

      { Transform basis by swapping component ComCtr and species SpeCtr }

      TransFormBasis(SpeCtr,ComCtr);

      LogX[ComCtr] := Conc[SpeCtr];
      PrecipCount := PrecipCount - 1;    { Dissolution }
      Conc[NowSpe1] := 0;

      NowSpe1 := NowSpe1 - 1;
      NowComTot := NowComTot + 1;
      NowComDiss := NowComDiss - 1;

   end;                          {End of Procedure DIssolve }






{ RestoreArrays returns the species and component concentrations to their
'proper' places in the Conc, LogX and TotErr arrays following a
solution involving precipitated solids. It also re-initializes the
arrays A, LogKZero, Tot, etc., to prepare for the next solution. }
   Procedure RestoreArrays;

   Var spe, com : integer;            { counters }
      ReConc : SpecArray;            { holds unscrambled Conc data }
      ReLogX, ReTotErr : Vector;     { holds unscrambled LogX, TotErr data }

   begin                                         { begin RestoreArrays }

      for com := 1 to S.NComp do begin            { for all current components }
         spe := CompSwaps[com];
         if (spe > 0) then begin                  { if originally a component... }
            ReLogX[spe] := LogX[com];
            ReTotErr[spe] := TotErr[com];
         end
         else ReConc[-1 * spe] := exp(LogX[com]) - TotErr[com]; { if precipitate }
      end;

      for spe := 1 to S.NSpec do begin            { for all current species }
         com := SpecSwaps[spe];
         if (com > 0) then ReConc[com] := Conc[spe]  { if originally a species }
         else begin                                { if a transformed component }
            com := com * -1;
            ReLogX[com] := ln(Conc[spe]);
            ReTotErr[com] := 0;
         end;
      end;

      Conc := ReConc;
      LogX := ReLogX;
      TotErr := ReTotErr;
      InitArrays(False);                        { Initialize arrays except LogX }

   end;                                        { end of RestoreArrays }




{**************************************************************************}
{ Now we begin the main algorithm for SolidSolve                           }
{**************************************************************************}

Begin                             { This is the main code for SolidSolve }

   SolidSolve := False;

   {If there are no solids which might precipitate, then Exit }
   If (S.STypeNum[Solid] = 0) then begin
      SolidSolve := Solve(S.CTypeNum[CTotDiss], S.STypeNum[Diss]);
      Exit
   end;

   SolidTries := 0;                 { number of tries at solution = 0 }
   PrecipCount := 0;                { number of solids precipitated = 0 }

   FillChar(CompSwaps,SizeOf(CompSwaps),0);     { initialize comp position }
   for i := 1 to S.NComp do
      CompSwaps[i] := i;

   FillChar(SpecSwaps,SizeOf(SpecSwaps),0);     { initialize spec position }
   for i := 1 to S.NSpec do
      SpecSwaps[i] := i;

   FillChar(Conc,SizeOf(Conc),0);

   NowSpe1 := S.STypeNum[Diss] ;      { number species in mass calculations }
   NowComTot := S.CTypeNum[CTotDiss];   { number of mass-limited components }
   NowComDiss := S.CTypeNum[CFixDiss];   { number of fixed activity components }

   Repeat                          { The case with possible diss. / precip. }

      Solve(NowComTot,NowSpe1);                        { Solve the system }
      If (not MassOK) then Exit;

      {ComCtr := DissolutionCheck;  }      { Check dissolution of solids }
      ComCtr := 0;                         {Dissolution of components impossible}
      {Now that I think about it, I do not }
      {understand why the dissolution routine}
      {ought to work or be useful }
      if (ComCtr <> 0) then Dissolve(ComCtr)
      else begin
         SpeCtr := PrecipitationCheck;  { Check precipitation of solids }
         if (SpeCtr <> 0) then Precipitate(SpeCtr);
      end;

      SolidTries := SolidTries + 1;              { Increment iterative counter }

   Until ((SolidTries=MaxCount) or ((SpeCtr=0) and (ComCtr=0)) or PhaseViol);

   If (PhaseViol or (SolidTries = MaxCount) or TooMany) then Exit;

   SolidSolve := True;

   If (SolidTries > 1) then RestoreArrays;      { unscramble arrays }

end;





{********************************************************************}
{   SweepSolve systematicaly varies a single component concentration,     }
{   (stored in the variable Xdata) solving the equilibrium for each value.}
{   It returns True if all problems are solved correctly, False if any    }
{   calculation fails to converge                                         }
{********************************************************************}

function SweepSolve: Boolean;

Var     ttr : integer;                       { Sweep point counter}
   com, spe : integer;                  {component and species counters}

{ Main body of procedure SweepSolve}

Begin
   SweepSolve := False;           {not solved yet}

   InitArrays(True);              { initialize vectors and arrays from S }
   ttr := 0;                      { set counter to zero }

   repeat                          { Main Sweep Loop }

      ttr := ttr + 1;                { increment point counter  }

      With S do
         case TiType of                { change the X variable in the Data Set }

           TiTotal: Comps[TitrX].TotalX := XData[ttr];  { Class 1 component sweep }

           TiFree: Comps[TitrX].GuessX := XData[ttr];  { Class 2 component sweep }

         end;                          { end of cases for change X }
      { end of With S }

      InitArrays(False);             {Now initialize all computational arrays, }
      {retaining the previous FreeX concentrations }
      {for all Known Total components }

      If (not SolidSolve) then Exit;  {if a calculation does not converge, Exit }
      {and return False}

      Replace;
      for com := 1 to S.NComp do     { record molarity of monitored comps, specs}
         YData[com,ttr] := S.Comps[com].Conc;
      for spe := 1 to S.NSpec do
         Ydata[spe + S.NComp, ttr] := S.Specs[spe].Conc;

   until (ttr = TitrPts) ;          { end titration loop }
   { since all calculations were successful, }
   SweepSolve := True;          { Return a value of True}

end;         {end of SweepSolve}

{********************************************************************}
{   TitrateSolve simulates a titration by solving the equilibrium problem     }
{   for a series of added titrant volumes (stored in the variable AddedVol).  }
{   It returns True if all problems are solved correctly, False if any        }
{   calculation fails to converge                                             }
{********************************************************************}

function TitrateSolve : Boolean;
Var    ttr : integer;                    {Titration point counter}
   com, spe : integer;               {component and species counters}
   VTotal : double;                  {Total volume for a given point}
   Dilution : double;                {Dilution factor during titration}
   Amount  : double;                 {Amount of titrated component, moles}

{ Main body of procedure TitrateSolve}

Begin
   TitrateSolve := False;          {not solved yet}

   InitArrays(True);               { initialize vectors and arrays from S }
   ttr := 1;                       { set counter to first point}

   repeat                          { Main Simulation Loop }

      VTotal := V0 + VolAdded[ttr];  {Calculate total volume of solution}
      Dilution := V0 / VTotal;       {Calculate the dilution factor}

      For com := 1 to S.NComp do     {Now for every component...}
      begin                        {Check to see if it has a total constraint}
         If ((S.comps[com].CType = CTotDiss)) then
         begin                   {if it does,...}
            If (com = TitrCation) then      {if this is present in the titrant}
            begin
               {Calculate the new total concentration, included added cation}
               Amount := OriginalSystem.Comps[com].TotalX * V0 +
               VolAdded[ttr] * CationConc;
               S.comps[com].TotalX := Amount / VTotal;
            end                  {End if com = Cation }
            else if (com = TitrAnion) then
            begin
               {Calculate the new total concentration, included added anion}
               Amount := OriginalSystem.Comps[com].TotalX * V0 +
               VolAdded[ttr] * AnionConc;
               S.comps[com].TotalX := Amount / VTotal;
            end
            else
               {Calculate the new total concentration by simple dilution}
               S.comps[com].TotalX := OriginalSystem.comps[com].TotalX * Dilution;
         end;                     {End if S.comps....}
      end;                       {End com for loop}

      InitArrays(False);             {Now initialize all computational arrays, }
      {retaining the previous FreeX concentrations }
      {for all Known Total components }

      If (not SolidSolve) then Exit;  {if a calculation does not converge, Exit }
      {and return False}

      Replace;

      for com := 1 to S.NComp do     { record molarity of monitored comps, specs}
         YData[com,ttr] := S.Comps[com].Conc;
      for spe := 1 to S.NSpec do
         Ydata[spe + S.NComp, ttr] := S.Specs[spe].Conc;

      ttr := ttr + 1;                { increment point counter  }

   until (ttr > SimuPts) ;          { end titration loop }
   { since all calculations were successful, }
   TitrateSolve := True;          { Return a value of True}

end;         {end of TitrateSolve}



                {End of Unit Solver}
end.
